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Path integral Monte Carlo with importance sampling for excitons interacting with an 
arbitrary phonon bath 
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The reduced density matrix of excitons coupled to a phonon bath at a finite temper- 
ature is studied using the path integral Monte Carlo method. Appropriate choices 
^ of estimators and importance sampling schemes are crucial to the performance of 

j^ the Monte Carlo simulation. We show that by choosing the population-normalized 

1^ estimator for the reduced density matrix, an efficient and physically-meaningful sam- 

^ pling function can be obtained. In addition, the nonadiabatic phonon probability 

I— I density is obtained as a byproduct during the sampling procedure. For importance 

Qh sampling, we adopted the Metropolis- adjusted Langevin algorithm. The analytic ex- 

?H pression for the gradient of the target probability density function associated with 

^ the population-normalized estimator cannot be obtained in closed form without a 

matrix power series. An approximated gradient that can be efficiently calculated 
^ is explored to achieve better computational scaling and efficiency. Application to a 

I~--~ simple one-dimensional model system from the previous literature confirms the cor- 

\0 rectness of the method developed in this manuscript. The displaced harmonic model 

Q system within the single exciton manifold shows the numerically exact temperature 

,__! dependence of the coherence and population of the excitonic system. The sampling 

_^ scheme can be applied to an arbitrary anharmonic environment, such as multichro- 

mophoric systems embedded in the protein complex. The result of this study is 
expected to stimulate further development of real time propagation methods that 
satisfy the detailed balance condition for exciton populations. 
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I. INTRODUCTION 

Recent 2D non-linear spectroscopy experiments suggested the existence of long-lived 
quantum coherence during the electronic energy transfer process within the Fenna-Matthews- 
Olson complex of green sulfur bacteria, marine algae and plants even under physiological 
conditional"^. These results attracted a large amount of attention from theoretical physicists 
and chemists. The energy transfer process usually has been modeled as the dynamics of 
excitons coupled to a phonon bath in thermal equilibrium within the single exciton man- 
ifold. This approximation leads to the famous spin-Boson Hamiltonian. The solution of 
this type of Hamiltonian has been studied extensively. For example, by assuming a certain 
relative magnitude between the reorganization energy and coupling terms, one can obtain 
quantum master equations valid in specific regime^^^. Another approximation, the Haken- 
Strobl-Reineker model works in both the coherent and incoherent regimes, but incorrectly 
converges to the high temperature limit in the long time even at the low temperature ^^. 
More recently, numerically exact approaches which interpolate both limits have been in- 
vestigated and applied to many systems of interest. Two of the most popular methods are 
the hierarchical equation of motion ISHlll a^^j ^j^e quasiadiabatic path integral methocP^'^. 
These methods are being actively developed, improved, and applied to many systems of 
interests^. 

Although having been successful in many applications, many of the models described 
above have assumed the phonon bath to be a set of independent harmonic oscillators and 
encode all the complexity of the bath environment in the spectral density, which is essentially 
a frequency dependent distribution of exciton-phonon coupling. However, for studying the 
anharmonic effects of a very sophisticated bath environment, like the protein complexes of 
photosynthesis, being able to directly include the atomistic details of the bath structure into 
the exciton dynamics has a distinct advantage. In other words, approaches that can evaluate 
the influence functional first suggested by Feynman and Vernorip^ have more straightforward 
descriptions and are applicable to arbitrary systems. Evaluation of the exact influence 
functional for arbitrary environment requires the simulation of the full quantum dynamics, 
which is still not practical with currently available computational resources. There have been 
several attempts to incorporate atomistic details of the large scale bath by combining the 
exciton dynamics and molecular dynamics simulations l^^"^. However, these theories are still 



in their early stages and the propagation scheme used does not satisfy some fundamental 
properties, like the detailed balance condition at finite temperature. In pursuit of more 
accurate theory, it is crucial to know the correct asymptotic behavior in the limit of infinite 
time. In this context, we decided to explore the numerically exact reduced density matrix 
in a finite temperature using path integral Monte CarlcP^'^ method. Recently, Moix et al 
applied path integral Monte Carlo for the equilibrium reduced density matrix of the FMO 
complex within the framework of open quantum systems^. 



II. THEORY 

A. Path integral formulation of the reduced thermal density matrix 

We want to evaluate the reduced density matrix of an excitonic system coupled to phonons 
on arbitrary Born-Oppenheimer surfaces at a finite temperature. For photosynthetic energy 
transfer, we usually restrict the excitons to be within the single exciton manifold because at 
normal light intensity, in average, one photon is present at a given time in the complexes of 
interest. However, the formulation itself is not limited to the single exciton manifold. The 
Hamiltonian operator for such a system can be written as 

H = ^ dR [Kn(-R) - Vg{R)] \m){m\ ® \R){R\ + ^ dR J„„(i?)|m)(n| ® \R){R\ 
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The Hamiltonian was written in terms of the diabatic basis \m^R) = \m) |-R), where m 
is the index for the exciton state and R is the phonon coordinate. V^(-R) is the potential 
energy surface (PES) of the phonons in the electronic ground state and Vm{R) is the PES 
of the phonons in the mth exciton state. T is the kinetic operator of the phonons defined as 
T = — yA^'^V^, where A^ is the mass tensor of the phonons. This expression is generally 
applicable to any molecular system with multiple potential energy surfaces. The reduced 
thermal density matrix ps is defined as the partial trace of the full thermal density matrix 



with respect to the bath degrees of freedom: 

P5 = ^TrBexp(-/3^; 

= ^ y rfi^o {Ro\ exp (-/3i7) |i^o), (2) 

where Z{P) is the partition function of the total system. We proceed by relying on the 
following identity: 

{Ro\exp{-/3H)\Ro) = {Ro\lexp (-j[] \ \Ro) 

dRi / dR2 • • • I uRm—i 



, I3H\ ,„,,„, / I3H\ ,„ 
X (iiolexp I -— I |i?M-i)(-RM-i|exp I -— I \Rm-2) ■ ■ ■ 

X {R^lexp i-j[] \Ri){Ri\exp (-J^j \Ro)- (3) 

For any positive integer M, the expression above is exact. When the Trotter decomposi- 
tion is applied, an imaginary timestep t = j^ is usually defined for convenience. Then, the 
thermal density matrix can be interpreted as an imaginary time evolution. In the limit of an 
infinitesimal imaginary timestep, the Trotter decomposition converges to the exact result, 

(iiilexp (-^J l^o) = (-Rilexp {-rH/h) \R^) 

= f dR2 [ dRs {Ri\e-^^^-^'^^\R^) 
X (i?3|e-"^^/''|i?2)(i?2|e-"'^<=''^/'''ii?o) + 0{t^). (4) 

Subsequently, we will recast the system part of i/cxc as a single matrix to simplify the 
notation, 

^exc = ^ dR E„,n{R)\m)in\ ® |i^)(i^|, 

, V^iR) - VJR) for m = n, 
Jmn{R) for m^n. 



With the single exciton manifold assumption, Emm corresponds to the optical gap of the 
m-th site. Now, the three terms in the integrand of the Eq. |4]can be written without Dirac 
notation, 

{R2\e'^''-'^''\RQ) = 6{R2 - Ro)e-^'"-'^'^^^^\ (6) 

where A = ^^. By the Eq. Q and Eq. |6 , 
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^-TEiRi)/2h^-TE(Ro)/2h _^ 0(t^). (7) 



Note that Eq. [7] is a matrix with the same dimension as the reduced density matrix of the 
system. Substituting Eq. [7] to Eq. |2| we obtain 

PS = y( n\ / '^-'^0 / dRi ■ ■ ■ cIRm-i 

^ ^-TE{Ro)l2ri^-TE{RM-i)/ri . . . ^-TEiRi)/h^-TE{Ro)/2h 

^ ^-TVg{Ro)/h^-rVg{Ri)/h . . . ^-rVg{RM^i)/ri 

^ g-{i?0--RA/-l)^A-l(fio-flAf-l)/4rg-(KM_i-fiM-2)^A-l(flM_i-fiM_2)/4r 



ujRq I dR\ ■ ■ • j cIRm—i 
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X ^-TE{Ro)/2h^-TE{RM-i)/h . . . ^-TE(Ri)/h^-TE{Ro)/2h 



PFIMC{R0,--- ,Rm-i) 

X _^-l3VpiMciRt),Ri,- ,Rm-i) fg\ 

K 

■^ V ' 

fgiRo,-,RM-i) 



where, 

^PIMC(^05 Rli ■ ■ ■ 5 Rm-i) = TT 2^ ^s(^i) 

j=0 

^-^ M 



+ 2_^ 'n~Q2k2^^^ ~ -'^mod{i+l,A/)} A^{-Ri " -Rmod(i+l,M)}- (9) 

i=0 ^ 



The expressions above show that the reduced thermal density matrix ps can be evalu- 
ated as an expectation value of ppimc(-Ro) " " " ,Rm~i) where the joint probability density 
function of the M iV-dimensional random variables {Rq, ■ ■ ■ , -Ra/-i) is fg- This type of mul- 
tidimensional integral can be efficiently evaluated using Monte Carlo integration. Because 
fg{Ro, ■ ■ ■ , Rm-i) is invariant to cyclic permutation of the phonon coordinate, usually the 
averaged estimator ppiMC over the cyclic permutation is used in the actual Monte Carlo 
evaluation: 

M-l 
PpiMC ~ T7 Z-^ PPlMciRi, -Rmod(i+l,Af), " " " , -Rmod(J+Af-l,M))- (10) 

4=0 

B. Population-normalized estimator and importance sampling 

In the previous approach described in Eq. [8} the phonon coordinates are sampled ac- 
cording the electronic ground state PES. The estimator should converge to the target 
quantity in the long time limit, taking into account the discretization error. As long as 
fg{Ro,--- ,Rm-i) is positive definite everywhere in the phonon space, the sampling effi- 
ciency depends on the selection of the probability density. Obviously, the actual distribution 
of the phonon coordinate depends heavily on the excited state PES. Therefore, the Monte 
Carlo points coordinates sampled according to the reduced dynamics of the bath by taking 
the partial trace with respect to the exciton degrees of freedom, as explored in multiple sur- 
face path integral Monte Carlo approaches, are expected to give the better estimates. This 
choice of the probability density reweights the estimator in the following way: 

fi{Ro, ■ ■ ■ , Rm-i) = Trs [PpimcI-'^O; ■ ■, -Rm-i)] fgiRo, ■ ■ ■ , Rm-i)-, 

r, (-D D ^ _ PpiMc(-Ro^ • • • ^ Rm-i) .. ,,, 

In the expression above, we call pi{Ro, ■ ■ ■ , Rm-i) the population normalized estimator for 
the reduced density matrix because the sum of its populations is always constrained to be 1. 
The effective energy gap term of —4 log Trp pj^g (-Ro, ■ ■ ■ , Rm-i) was added to the Eq. p^to 
enable the phonons follow the excited state dynamics depending on the exciton state ps- For 
the estimator of the reduced density matrix in Eq. |8| the normalization must obtained by the 
estimates of its diagonal elements, leading to more uncertainties in the coherence. However, 
the population-normalized estimator preserves the correct normalization by construction, 
and does not suffer from any additional uncertainty. 



Local gradient information can improve the efficiency and scaling of the sampling pro- 
cedure by means of a gradient-based approach such as the Metropolis-adjusted Langevin 
algorithm (MALA).'^^'^ However, the exact closed form of the gradient of the effective en- 
ergy gap term, log TrsP pjMc (-^0; ' ' ' > Rm-i) can only be expressed as a function of a power 
series of matrices. Nevertheless, with the following approximation: 



n " 1 / \ 



(12) 

fc=0 k=0 

an accurate approximated of the gradient can be obtained and employed in the sampling 
procedure. 
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Here, Vj is the gradient operator with respect to Ri. 

Note that if we choose an appropriate Metropolis criterion, no bias in the distribution is 
introduced even with the approximate gradientpS. Firstly, a trial move R[ obtained by 

R[ = Ri + /ii(i?o, • • • , -Rm-i) At + ^iVAt, (14) 

where At is the timestep for the Monte Carlo step and ^i is a A^-dimensional vector of 
independent standard Gaussian random variables. Then, R'^ is probabilistically accepted 
according to the acceptance ratio. 
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The Monte Carlo timestep At is only a tunable parameter for the Monte Carlo sampling 
procedure and not related to the physics of the simulated system. 



Parameters 


Value 


fell 


4 X 10-5 


^22 


3.2 X 10-5 


j;ii 


7 


2;22 


10.5 


en 





£22 


2.2782 X 10-5 


c 


5 X 10-5 


a 


0.4 


3^12 


8.75 


m 


3.6743 X 103 



TABLE I. Summary of the parameters for the model system by Alexander et aP^. All values are 
given in atomic units. 

III. APPLICATION 

A. Alexander's ID test model 

Our formulation is equivalent to the multiple electronic state extension of matrix multi- 
plication path integral (MMPl) method of Alexander^^'^ when the population normalized 
estimator is chosen and only the vibrational degrees of freedom are considered. Therefore, 
the ID model employed in Ref. [30] was calculated to test the validity of our method. The 
elements of the electronic Hamiltonian in this model are given by, 

^22(3;) = -k22{x - X22f + £:22, 

^12(3^) = cexp [— a(x — Xi2)^] , (16) 

The total nuclear probability density evaluated as histograms from the Metropolis random 
walk and MALA simulations are compared to the grid-based result from Alexander et al^ 
in Fig. [ij The distributions converged to the exact probability density after 2 x 10^ steps 
with 8 beads at both temperatures of 8K and 30K. 
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FIG. 1. The estimated nuclear probability densities of Alexander's modeP^lat (a) 8K and (b) 30K. 
For path integral Monte Carlo simulations, densities were obtained by histograms with 50 bins. 
The discretization number of 8 was enough to converge to the exact probability densitiies. 



B. Model of a chromophore heterodimer with displaced harmonic oscillators 



To test the proposed method, a system of two chromophores in a photosynthetic complex 
was modeled using displaced harmonic oscillator model. In this model, the ground and 
excited electronic states of the monomer are modeled as harmonic oscillators with different 
displacement, but the same harmonic constantP. The thermal reduced density matrix was 
calculated within the single exciton manifold. The Hamiltonian for this model is then given 



Parameter 


Value 


ki 


2.227817 X 10-3 


k2 


2.227817 X 10-3 


di 


3.00000 


d2 


2.00000 


ei 


8.064745 X 10"^ 


£2 


7.976238 x lO'^ 


J 


-4.738588 x 10"^ 


7711 


3.418218 X 10^ 


1712 


3.418218 X 10*^ 



TABLE II. Summary of the parameters for the displaced harmonic oscillator model used in 



Sec. IIIB All values are given in atomic units. 



as follows: 

1 



Vg{xi,X2) = ^{kixl + k2xl), 



Ve(Xi,X2) 



^ki{{xi - di)^ - xf} + El J 

J Ik2{{x2 - d2y - xl + 62} 




M={ ' . (17) 



Some of the parameters were set according to our molecular dynamics/quantum chemistry 
calculation of the FMO complex^. The parameter values are listed in table [nj 

The model system was simulated at seven different temperatures ranging from 30K to 
300K with a number of beads (discretization number) of 4, 8, 16, 32 and 64. The number 
of timesteps propagated in each simulation was 4 x 10^. The value of each timestep was 
tuned so that the acceptance ratio of the MALA run is close to 0.574, and 0.234 for the 
Metropolis random walk as maintaining these acceptance ratio is known to provide most 
efficient sampling^. We used non-overlapping batch means^ with a batch size of 10® to 
estimate the standard error of the correlated samples. The batch size was adjusted so that 
the null hypothesis of uncorrelated batches was not rejected by using Ljung-Box test^- at a 
significance level of 5%. 
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(c) MALA 

Q = 21.2888 
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(d) Random Walk 
Q = 17.4696 
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FIG. 2. Estimates of (1,2) matrix elements of the thermal reduced density matrix evaluated 
using MALA and Metropolis random walk at 77K with 64 beads. MALA estimate has a smaller 
confidence interval thus a more accurate estimate than that of the Metropolis random walk. The 
error bar indicates the 95% confidence interval evaluated with the batch means. The 0.95 quantile 
of the x^ distribution with 13 degrees of freedom is 22.362 and both Ljung-Box statistics (Q) are 
smaller. Thus, the uncorrelation hypothesis is not rejected in both cases at the 5% significance 
level. 

As shown in Fig. |2| the standard error of the simulation decreases modestly as the number 
of Monte Carlo steps increases. Fig. [3] shows the temperature dependence of the estimates 
of reduced density matrix elements as a function of various discretization numbers using 
MALA. Although the Metropolis random walk simulation gives a smaller confidence interval 
for the 4 bead case, MALA provides better estimates as the dimension of the sample space 
increases. The Metropolis random walk result is given in Fig. |4} While the population of 
the low energy site decreases as the temperature increases, the quantum coherence does not 
monotonically decrease. We believe that this pheonomenon is an artifact of an insufficient 
discretization number at low temperatures. As can be seen in Fig. |3| 64 or more beads 
are needed for the coherence to converge at 77K, while 16 beads are enough at 300K with 
acceptable accuracy. This is a well known limition of imaginary time path integral Monte 
Carlo simulations. Figure [5] shows the probability density function of the phonon coordinate 
at 77K and 300K. The population difference in the reduced density matrix is reflected to 
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FIG. 3. Estimates of matrix elements of the thermal reduced density matrix evaluated at 30K, 
50K, 77K, 140K, 225K and 300K with different discretization numbers of 4, 8 and 16 using MALA. 
(a) is the (1,1) element, (b), (c) and (d) are (1,2), (2,1) and (2,2) elements, respectively. The error 
bar indicates the 95% confidence interval evaluated with the batch means. 

the difference in the probability mass of the two diabatic potential energy minimum at (3, 0) 
and (0,2). 

IV. CONCLUSION 

We explore a method for obtaining the thermal reduced density matrix of an exciton 
system coupled to an arbitrary phonon bath for path integral Monte Carlo simulation. 
Note that our scheme is closely related to the path integral Monte Carlo simulation for 
nonadiabatic systems for vibrational coherenc d^°*^^ * ^^ l Although the phonon state can be 
obtained as a byproduct, we mainly focused on the evaluation of the reduced density matrix 
of the excitonic system to explore the asymptotic behavior of the populations and coherences 
in this paper. In addition, we implemented an importance sampling scheme for better spatial 
scaling and sampling efficiency. Although the path integral Monte Carlo cannot evaluate the 
real time evolution of density matrices, the method gives the exact asymptotic values with all 
quantum effects from both the system and bath environments if a sufficient number of beads 
are used. We believe that in some of the cases where the bath has a nontrivial coupling 
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FIG. 4. Estimates of matrix elements of the thermal reduced density matrix evaluated at 30K, 
50K, 77K, 140K, 225K and 300K with different discretization numbers of 4, 8 and 16 using random 
walk Metropolis, (a) is the (1,1) element, (b), (c) and (d) are (1,2), (2,1) and (2,2) elements, 
respectively. The error bar indicates the 95% confidence interval evaluated with the batch means. 
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FIG. 5. The phonon probability density function evaluated at (a) 77K and (b) 300K with 16 beads 
using MALA. At the lower temperature, the contribution of the exciton with lower energy at (0, 2) 
becomes larger. Therefore, the population differenece becomes more distinct, as can be seen in the 
temperature dependence of the exciton population in Fig. [3j 
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to the system, or the non-Markovianity of the bath manifests very strongly, treating the 
environment around the system of interest as a set of harmonic oscillators is not sufficient. 
If this is the case, the system should be studied in its entirety. We are trying to develop a 
real time propagation scheme to treat the system exactly, and the bath semiclassically. The 
method studied in this paper offers a foundation for it by providing the correct asymptotic 
behaviors. 
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